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The flow behaviors of polymer melt composed of short chains with ten beads between parallel 
plates are simulated by using a hybrid method of molecular dynamics and computational fluid 
dynamics. Three problems are solved: creep motion under a constant shear stress and its recovery 
motion after removing the stress, pressure-driven flows, and the flows in rapidly oscillating plates. In 
the creep/recovery problem, the delayed elastic deformation in the creep motion and evident elastic 
behavior in the recovery motion are demonstrated. The velocity profiles of the melt in pressure- 
driven flows are quite different from those of Newtonian fiuid due to shear thinning. Velocity 
gradients of the melt become steeper near the plates and flatter at the middle between the plates 
as the pressure gradient increases and the temperature decreases. In the rapidly oscillating plates, 
the viscous boundary layer of the melt is much thinner than that of Newtonian fluid due to the 
shear thinning of the melt. Three different rheological regimes, i.e., the viscous fluid, visco-elastic 
liquid, and visco-elastic solid regimes, form over the oscillating plate according to the local Deborah 
numbers. The melt behaves as a viscous fluid in a region for wr^ < 1, and the crossover between 
the liquid-like and solid-like regime takes place around ujt" ~ 1 (where ui is the angular frequency 
of the plate and and are Rouse and a relaxation time, respectively). 
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I. INTRODUCTION 

Many products in our daily life contain soft matters 
(e.g., polymeric liquids, colloidal dispersion, and liquid 
crystals). Soft matters involve complex internal degrees 
of freedom, such as orientation of molecules, dispersed 
particles, and phase-separated structure, and they ex- 
hibit peculiar flow behaviors coupled to the micro-scale 
dynamics, e.g., visco-elastic flow, shear thinning or thick- 
ening, and flow-induced phase transition. Investigation 
of the complicated flow behaviors of soft matters are 
of great importance in various science and engineering 
flelds, such as fluid mechanics, materials science, biolog- 
ical science, mechanical engineering, and chemical engi- 
neering. To predict the flow behaviors of soft matters 
by computer simulation is challenging from both an aca- 
demic and a practical point of view. The difficulty in 
simulating soft-matter flows also is also due to the com- 
plicated couplings between the microscopic dynamics of 
internal degrees of freedom and macroscopic flow behav- 
iors. In the present paper, we demonstrate that multi- 
scale modeling is a promising candidate for soft-matter 
simulations. We solve the complicated visco-elastic flow 
behaviors of polymer melt by using a multiscale hybrid 
simulation method. 

Usually, when computer simulations of polymer melt 
flows are performed, either computational fluid dynam- 
ics (CFD) or molecular dynamics (MD) are employed. In 
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the case of CFD, mechanical properties of fluids must be 
modeled mathematically in advance as a form of "con- 
stitutive relationship" to be used in simulations. CFD is 
thus valid only for the cases in which mechanical proper- 
ties of fluids are not too complex. Polymer melts, how- 
ever, have very complicated mechanical properties in gen- 
eral. In the case of MD simulation, in contrast, fluids con- 
sist of huge numbers of molecules of arbitrary shapes. In 
principle, MD simulation is thus applicable for any flows 
of any complex fluids. However, the drawback here is 
that enormous computational time is required to resolve 
the dynamics of all the constituent molecules. Hence, 
MD simulation is not yet applicable to problems which 
concern large-scale motions far beyond molecular size, as 
is done in the present paper. In order to overcome the 
weaknesses of the individual methods, we have developed 
a multi-scale simulation method that combines MD and 
CFD. 

In our hybrid simulation method, the macroscopic flow 
behaviors are solved using a CFD scheme, but, instead of 
using any constitutive equations, a local stress is calcu- 
lated by using a non-equilibrium MD simulation associ- 
ated with each lattice node of the CFD simulation. The 
basic idea of this type of hybrid simulation method was 
flrst proposed by E and Engquist^iS^ where the heteroge- 
neous multiscale method (HMM) is presented as a gen- 
eral methodology for the efficient numerical computation 
of problems with multiscale problems. The HMM has 
also been applied to the simulations of complex fluidsi^ 
The equation-free multiscale computation proposed by 
Kevrekidis et al. is also based on a similar idea and has 
been applied to various problemsi^ De et al. have de- 
veloped a hybrid method, which is called scale bridging 
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method in their paper, that can correctly reproduce the 
memory effect of a polymeric liquid and demonstrate the 
non-linear visco-elastic behaviors of polymeric liquid be- 
tween oscillating plates. - The methodology of the present 
hybrid simulation is the same as that of the scale bridg- 
ing method. The present hybrid method is also similar 
to the CONNFFESIT approacbJ'^'^ in that the macro- 
scopic local stresses are obtained from the configurations 
of polymer molecules. In the present method, however, 
we couple a CFD simulation with an MD simulation 
which involves full interactions of each molecule rather 
than the stochastic dynamics of single model polymers 
lacking any direct interactions between molecules. Thus 
the present method can reproduce flow behaviors caused 
by the many-body effect, which is quite important, es- 
pecially in the glassy polymers. In previous work^'i^, 
we investigated the efficiency and accuracy of the hybrid 
method. We also clarified the property of noise arising 
in the hybrid method by comparing the results of the 
hybrid method with those from fluctuating hydrodynam- 
ics. In addition, we analyzed the complicated rheological 
properties of a supercooled polymer melt in the viscous 
diffusion layer arising over the rapidly oscillating plate 
by using the scale bridging method. 

In the present paper, various flow behaviors of poly- 
mer melts between the parallel plates are simulated using 
the scale bridging method. Creep motion under a con- 
stant shear stress and its recovery motion in the stress- 
free state, pressure-driven flows, and flows in rapidly os- 
cillating plates are investigated. In the creep/recovery 
simulation, the evident elastic motion of polymer melt 
is demonstrated. The result is also compared with that 
from a model constitutive relation. For pressure-driven 
flow, the velocity profiles of melts with various temper- 
atures and pressure gradients are demonstrated. The 
shear-thinning behaviors of melts are also investigated. 
In the problem of oscillating plates, the viscous bound- 
ary layer arising over the oscillating plate and rheologi- 
cal properties of the boundary layer are investigated for 
various temperatures, frequencies of the plate, and am- 
plitudes of strain of the system. The parameter sets used 
in the present paper are different from those used in the 
previous paper^^; in the previous paper, the amplitude 
of velocity of the plate is specified as a parameter instead 
of the amplitude of strain of the system in the present 
paper. 

For the problems of pressure-driven flow and oscillat- 
ing plates, the velocity profiles are quite non-uniform be- 
tween the plates, and two different characteristic length 
scales appear in each problem that must be resolved: one 
is that of a polymer chain, which is the scale of the MD 
simulation, and the other is that of the flow behaviors 
of the melt, e.g., the width between the plates or the 
thickness of boundary layer arising over the oscillating 
plate, which is the scale of the CFD simulation. These 
problems constitute important applications of multiscale 
modeling; it is quite difficult to solve each problem by us- 
ing a full MD simulation because the width of the plates 



and thickness of the boundary layer are much larger than 
the size of a polymer chain. 

In the following text, the multiscale modeling and sim- 
ulation method for the present problem are stated in Sec. 
II. Some model constitutive relations for polymeric liq- 
uids and the difficulties with using these in CFD simu- 
lations are also mentioned in Sec. II. In Sec. Ill, the 
simulation results for the creep and recovery, pressure 
driven flows, and oscillating plates are given. The rheo- 
logical properties are also discussed in Sec. III. Finally, a 
summary and an outlook for the hybrid simulations are 
given in Sec. IV. 



II. MULTISCALE MODELING AND 
SIMULATION METHOD 

We consider the polymer melt with a uniform density 
Po and a uniform temperature Tq between two parallel 
plates (see Fig. [2a)). The upper- and lower- plate can 
move in the x-direction. The melt is composed of short 
chains with ten beads. The number of bead particles 
composing each chain in the MD simulation is repre- 
sented by iVb- Thus iVb = 10. AU of the bead parti- 
cles interact with a truncated Lennard-Jones potential 
defined byi^, 



4e [((7/r)i2 - (a/rf] + e {r < l^'^a), 
(r > 2^l^a). 



(1) 



By using the repulsive part of the Lennard-Jones poten- 
tial only, we may prevent spatial overlap of the particles. 
Consecutive beads on each chain are connected by an 
anharmonic spring potential. 
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fcei?^ln[l-(r/i?o)'] 
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where fcc=30e/(T^ and RQ — X.ha. The number density of 
the bead particles is fixed at pQ/m=l/a^, where m is the 
mass of the bead particle. With this number density the 
configuration of bead particles becomes severely jammed 
at a low temperature, resulting in the complicated non- 
Newtonian viscosity and long-time relaxation phenomena 
typically seen in glassy polymersiJ^jii 

We assume that macroscopic quantities are uniform in 
the X- and z-directions, d/dx^d/dz=0. The macroscopic 
velocity Va is described by the following equations, 



dt 



dy 



(3) 



and Vy=Vz=0, where t is the time and Gap is the stress 
tensor. For pressure-driven flow, — AP is also added to 
the right-hand side of Eq. ([3]), where AP is a uniform 
pressure gradient in x-direction. Here and afterwards, 
the subscripts a, and 7 represent the index in Carte- 
sian coordinates, i.e., {a^P ,^}={x,y . We also assume 
the non-slip boundary condition on each plate. 
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FIG. 1: Schematics for geometry of problem, mesh system, and time-evolution scheme 
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The constitutive relation of the stress tensor may be 
written in a functional of the history of velocity gradients 
on the fluid element, 

(Jc^pit^xo:) ^ F^p[na,p{t' ,x'^{t')% with t' <t, (4) 

where K^p is the velocity gradient, Kap = dva/dxp, and 
x'^{t') represents the path line along which a fluid element 
has been moving. Note that the temporal value of the 
stress tensor of a fluid element depends upon the previous 
values of the velocity gradients of the fluid element. In 
the one-dimensional problem, however, the equation sim- 
plifies to a functional involving only the local strain rate 
i—dvx/dy since the macroscopic velocity is restricted to 
the x-direction where the macroscopic quantities are as- 
sumed to be uniform; 



=F,y[7(i',y)], with t'<t. 



(5) 



Constructing a model for the constitutive relation is one 
of the most important goals in complex fluid research, 
and many models have been proposed . ""^^'^^ The simplest 
model for the visco-elastic motion is one proposed by 
Maxwellii, which is written as 



a^y + A 



mi- 



(6) 



where Ai is a time constant and rji a viscosity. For the 
steady state, this equation simplifies to the Newtonian 
fluid with viscosity ryi, while for sudden changes in stress, 
the time derivative term dominates the left-hand side of 
the equation, and then integration with respect to time 
gives the Hookean solid with elastic modulus G = 771 /Ai. 
The model which includes the time derivative of 7 in the 
right-hand side of Eq. ^ has also been proposed by 
JeffreySfi^ 



f^xv + A 



daxy 



= i]i 7 + A: 



dj 
■~dt 



(7) 



This equation contains two time constants Ai (the "relax- 
ation time" ) and A2 (the "retardation time"). This model 
may reproduce the delayed elastic motion for the sudden 
change of stress. Models also exist for the multiple relax- 
ation modes and nonlinear simulations involving either 
the shear-rate dependent viscosity 771(7) or the nonlinear 



stress term/iSiSS Usually, when one performs the CFD 
simulations for the polymer melt flows, one chooses an 
appropriate model of them empirically according to the 
physical properties of the problem or the computational 
convenience. The time constants, viscosities, and their 
dependence on the shear rate in the equation must be 
specified in advance. However, no systematic methods 
exist to choose a suitable model for each problem, and 
determination of time constants, viscosities, and their de- 
pendence on the shear rate is quite difficult, especially for 
glassy materials. In our multiscale modeling, instead of 
using explicit formulas for the constitutive relation, the 
local stress is generated by the non-equilibrium MD sim- 
ulation associated with each local point. See Fig. [T] (b). 

We use a common finite volume method with a stag- 
gered arrangement for the CFD calculation, where the 
velocity is computed at the node of each slit and the 
stress is computed at the center of each slit. For the 
time-integration scheme, we use the simple explicit Eu- 
ler method with a small time-step size At. The local 
stress is calculated at each time step of CFD by using 
the non-equilibrium MD simulation with a small cubic 
MD cell associated with each slit according to a local 
strain rate. In each MD simulation, we solve the so- 
called SLLOD equations of motion with the Gaussian 
iso-kinetic thermostat:^^-^^ 



— ^ = ^^+7i?4e, 
at m 



(8a) 



(8b) 



where is the unit vector in the x direction and the 
indexes k and j represent the fcth polymer chain (fc = 
1, • • ■ , A^p) and the jth bead (j = 1, • • • , A^b) on each 
chain, respectively. Here the number of polymer chains 
contained in each MD cell is represented by Np. Rj and 

p'^ +m'jRyjex are the position and momentum of the jth 

bead on the fcth polymer chain, respectively, fj is the 
force acting on the jth bead on the kth polymer due to 
the potentials described in Eqs. ([T]) and and 7 is the 
shear rate subjected on each MD cell, which corresponds 
to the local shear rate at each slit of the CFD calculation. 
Note that, in the SLLOD equations, /m represents the 
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deviation of velocity of each particle from the mean flow 
velocity jRyjE^ in the MD cell. The friction coefhcient 
C is determined to satisfy the constant temperature con- 
dition dT/dt = with T = Y,j.k IPj iVSmTVpTVb, where 
fc represents the summation all over the bead particles 
in each MD cell. The friction coefficient ( is calculated 
as 

j,k j.k 

We integrate Eq. ^ with Eq. ^ by using the leapfrog 
algorithmi^ The space integral of the microscopic stress 
tensor reads as 

k—1 j — 1 allpairs 



EE 

k=ij=i 



dC 



(10) 



where we rewrite the momentum of the jth bead on the 
fcth chain, p'^ + m^Ry^.ex, as p^. ^ in the right-hand 



side of Eq. (?? represents the relative vector R 



Rf 



between the two beads, R^ and R^ , in the second term 
and the relative vector R^ — -R^+i between the two con- 



secutive beads on the same chain, Rj and R'j+i, in the 
third term. 

In the present problem we cannot assume a local equi- 
librium state at each time step of the CFD simulation 
since the relaxation time of the stress may become much 
longer than the time-step size of the CFD simulation (in 
which the macroscopic motions of the system should be 
resolved). In the current simulations, the simple time- 
averages of the temporal stresses of the MD (averaged 
over the duration of a time-step of the CFD simulation) 
are used as the stresses at each time step of the CFD cal- 
culation without ignoring the transient time necessary for 
the MD system to be in steady state. See Fig. [Tfc). Thus 
the time integration of the macroscopic local stresses aa/3 
are calculated with the microscopic stress tensor Eq. 
as 

/.t+At 

^a/3{t, y) = J cTafsit', y)dt' 



/ n^p{T;^{t,y))dT, (11) 



"MD 

where Ai is the time-step size of the CFD simulation and 
Zmd is the side length of the cubic MD cell. Note that 
the second argument of IIq,^ in Eq. (jlip . i.e., 7(i, y), is 
constant in the integral interval. This indicates that the 
shear rate to which each MD cell is subjected is also con- 
stant over a duration Ai in each MD simulation. The 
final configuration of molecules obtained at each MD cell 
is retained as the initial configuration for the MD cell at 
the next time step of the CFD. Thus we trace all tem- 
poral evolutions of the microscopic configurations with 



a microscopic time step so that the memory effects can 
be reproduced correctly. Note that compared to a full 
MD simulation, we can save computation time with re- 
gard to the spatial integration by using MD cells that 
are smaller than the slit size used in the CFD simulation. 
The efficiency of the performance of the present hybrid 
simulation is represented by a saving factor defined by 
the ratio of the slit size used in the CFD simulation Ax 
to the cell size of the MD simulation /md, Acc/Zmd- It 
also should be noted that, in addition to the saving factor 
Ax/Zmd, the present hybrid method is quite suitable as 
a parallel computational algorithm since the MD simula- 
tions associated with each mesh of the CFD, which cost 
a large part of the total simulation time, are performed 
independently. 

We solve the polymer melt flows in the following three 
problems by using the hybrid method; (i) creep motion 
under a constant shear stress and recovery after remov- 
ing the stress, (ii) pressure-driven flows and (iii) flows 
in rapidly oscillating plates. In the first problem, the 
simple visco-elastic behavior of the melt is demonstrated 
in order to verify that the present hybrid method can 
reproduce the visco-elastic motion correctly. In the sec- 
ond problem, we demonstrate the peculiar velocity pro- 
file of the melt produced due to shear thinning near the 
plate. We also investigate the local rheology and the 
microscopic configurations at local points. In the third 
problem, the viscous boundary layers arising near the os- 
cillating plates are resolved. The macroscopic quantities 
are quite non-uniform in the boundary layers, and the 
local rheological properties are also changed drastically 
according to the local flow fields. 



III. SIMULATION RESULTS 



Hereafter, the quantities normalized by the units of 
length (7 and time tq=^ ma'^ j t are denoted with a hat 
" " " . In the following simulations, we fix the time-step 
size of the CFD simulation At, sampling duration of the 
MD simulation tMD and time-step size of the MD sim- 
ulation Ar as At=£MD = l and Af=0.001, respectively. 
Thus, 1000 MD steps (M = 1000 in Fig. ^c)) are per- 
formed in each MD cell at each time step of the CFD 
computation. One hundred chains with ten beads are 
confined in each cubic MD cell with a side length ^md=10; 
thus iVp=100 and 7Vb=10. The density of the melt is 
fixed to p—\. The temperature of the melt is Tb=0.2 
in the creep and recovery problem, and To=0.2 and 0.4 
are considered in the problems of pressure-driven flow 
and oscillating plates. At this number density and a low 
temperature, the polymer melt involves complicated non- 
Newtonian rheology. 
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FIG. 2: The time evolution of the strain of the system, F — [ux{y — 2H) ~ Ux(y = 0)] /2H. The solid line shows the creep 
motion of the melt. The dashed line shows the recovery motion of the melt after removing the stresses on the plates at f=20000. 
The dotted line shows the result for the Newtonian fluid. The dash-dotted line shows the result given by the Jeffreys model in 
Eq. (O. The inset shows the time evolution of the strain in the beginning of creep motion. 



A. Creep and Recovery 

The upper and lower plates were each subjected to 
a constant shear stress do at a time t—Q. For a pos- 
itive value of (To J the upper plate slides right and the 
lower plate left in Fig. [Tfa). We set the shear stress o-q 
and the distance between the plates 2H as o-o=0.1 and 
2.ff=800, respectively. The distance between the plates 
2H is divided into seventeen slits. The present mesh sys- 
tem is a little different from that in Fig. [TJb). In the 
present mesh system, sixteen MD cells are used to sam- 
ple the local stresses; these cells are associated with the 
nodes of the slits and the local velocities are calculated 
at the centers of each slit. Thus the saving factor in the 
present simulation is Ax/Zmd — 4.7. The stress to which 
the plate is subjected is transmitted inside the material 
immediately, and after a short transient time, the local 
macroscopic quantities of the melt are in uniform states 
between the plates (although they fluctuate due to the 
noise arising from each MD cell). The transient time 5t 
may be estimated from Fig. ^a,s 5t ^ 500. The system 
is deformed uniformly after a short transient time. 

We measure the strain of the system T{t) as 

P(^^^ ^ Uxiv = 2iJ,^) - u^{y ^ 0,t) 

where Ux{y,t) is the displacement of the melt in x- 
direction, which is calculated as 

Ux{y,t) = / v^{y,T)dT. (13) 
Jo 

Figure [2] shows the time evolution of the strain of the 
system T{t). In the flgure, the solid line shows the creep 
motion and the dashed line shows the recovery motion af- 
ter removing the stresses on the plates at a time ^=20000. 



The inset zooms in on the creep motion. It can be seen 
that, in the creep motion, the strain F increases rapidly at 
the beginning, say i ~ [0, 5000], and, as time passes, the 
time evolution of the strain becomes linear. The rapid 
evolution of the strain at the beginning of the creep mo- 
tion is caused by the elasticity of the meld. However, 
since it is different from the pure elastic deformation, 
in which a finite strain is instantaneously obtained when 
subjecting the stress, the strain of the melt evolves rather 
continuously. This behavior demonstrates the delayed 
elastic deformation of the melt. The linear evolution of 
the strain corresponds to the viscous flow. The viscosity 
of the melt in the linear evolution regime rji is calculated 
from the slope, which reads ?7i=1007. The result of the 
Newtonian fluid with the viscosity 7)1 is also shown in 
Fig [5] by the dotted line for comparison. When remov- 
ing the stress on the plate (the dashed line in Fig. [2]), 
the strain of the system starts to decrease very slowly to 
the line of the Newtonian fluid. This recovery motion of 
strain demonstrates the evident elastic motion; the elas- 
tic strain of the melt stored in the creep motion is being 
released after removing the stresses on the plates. 

The delayed elastic deformation may be described us- 
ing the Jeffreys model shown in Eq. ([7]). By integrating 
Eq. ([7]) with (Jxy — o'o8(t), where Q{t) is the step func- 
tion, the time evolution of the strain in the creep motion 
is written as follows: 

r(0 = ^(l-e-*)+^i, (14) 

where Gi = 7]i/(Ai — A2). The first term in the right-hand 
side of Eq. represents the delayed elastic deforma- 
tion and the second term the viscous flow. By fitting 
our result to Eq. (fT4|) . the viscosity ryi, elastic modulus 
Gi, and retardation time A2 are estimated as 7yi=1007, 
Gi=0.38, and A2 = 13003. The result of the Jeffreys model 
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TABLE I: Parameter values for pressure-driven flows and 
simulation times. T is the temperature and AP is the pressure 
gradient of the melt. tMax is the simulation time for each case. 
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FIG. 3: The time evolution of the standard deviation between 
the local stresses and the subjected stress on the plates ctq 
described in Eq. (|16|) . 



with these parameter values is shown in Fig. [2] (inset) as 
a dash-dotted line. It is evident that the creep motion ob- 
tained by the hybrid simulation can be fitted well by the 
delayed elastic motion described by the Jeffreys model 
Eq. p4|l . The recovery motion for the Jeffreys model is 
written as 

r(t)=(7o/Gi (1-6-^^)6""^+ Too, (15) 

where T^a is a elastic strain remaining at an infinite time. 
Eq. (|15p is obtained by integrating the strain rate r(t) 
(for t > t'), which is obtained by integrating Eq. ([7]) 
with axy = cro[l - Q{t - t')] using f (to) for Eq. (fT4| . 
from a time t to the infinity t = oo. In the recovery 
motion, however, the parameter values in the Jeffreys 
model differ from those in the creep motion since the 
mechanical property of the melt changes according to the 
state of motion. The result by the hybrid simulation in 
the recovery motion may be fitted to Eq. (jlSp with Gi ~ 
0.14, A2 = 1.6 X 10", t' = 20000, and Too^l-TG which 
is a value of strain obtained by the hybrid simulation 
at t = 2.5 X 10". The result of the Jeffreys model with 
these parameter values in recovery motion is shown as 
the dash-dotted line in Fig. [51 

As previously mentioned, the local macroscopic quan- 
tities fluctuate due to the noises arising from each MD 
cell. These fluctuations can be measured by the standard 
deviation of the local stresses from the subjected stress 
on the plates, which is defined by 



The time evolution of the standard deviation S with 
6=100 is shown in Fig. [31 After a short transient time 
5t, St ~ 500, the local stresses fluctuate around the uni- 
form state with a subjected stress on the plates ctq- The 
standard deviation in the fluctuating states is estimated 
to be approximately E ~ 0.045. 



B. Pressure-driven flow 

We consider a uniform pressure gradient in the x- 
direction AP. Both the upper- and lower-plates are at 
rest. For the present problem, Eq. ^ is modifled to 
include the uniform pressure gradient; — AP is added in 
the right-hand side of Eq. The distance between 

the plates is set as 2H = 1600. The mesh system is 
the same as that shown in Fig. [IJb). We use thirty-two 
slits in the distance 2H. The saving factor in this case is 
Ax/Imb = 5. Hybrid simulations are performed for the 
melts using various parameter values and arc shown in 
Table [H 

Figure [H shows the comparison of the velocity profiles 
normalized by the velocities at the middle between the 
plates Vb, Vq — Vx{y — H), for the melt and for the 
Newtonian fluid. The velocity proflles of the melt with 
T=0.4 are time-averaged over i = [50001,70000] and 
those of the melt with T=0.2 are time-averaged over t = 
[200001,400000] for AP=1 x 10"*, i = [100001,200000] 
for AP=3 X 10""*, and i = [50001, 100000] for AP=5 x 
10^"*. The value of Vq for each flow is also shown in Table 
im It can be seen that the velocity proflles of the melts 
are quite different from those of the Newtonian fluid. As 
the pressure gradient increases, the velocity gradients of 
the melts become steeper near the plates and are more 
gradual (or rather a plateau for T — 0.2) in the middle 
region. This feature is enhanced at a low temperature 
T = 0.2. The dependence of the velocity at the middle 
between the plates Vq on the pressure gradient AP for 
the melt is quite different from that for the Newtonian 
fluid. For the Newtonian fluid, the maximum velocity Vq 
depends linearly on the pressure gradient, while for the 



TABLE II: The velocity of the melt at the middle between 
the plates Vb in pressure-driven flows. For the Newtonian fluid 
with a viscosity i/q, Vq is written as Vb = —{APpo/2uo)H^ . 
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= 0.2 


f 


= 0.4 


AP X 


10* Vb 


AP X 


10* Vb 


1 


0.024 


3 


1.67 


3 


0.31 


5 


3.90 


5 


1.84 
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(a) f = 0.2 (b) f = 0.4 

FIG. 4: Comparison of the normalized velocity profiles in the steady states of the melt to those of the Newtonian fluid for the 
temperature T = 0.4 [for (a)] and T = 0.2 [for (b)]. Vb is the velocity at the middle between the plates. The value of Vb for 
each flow is shown in Table |ll] Note that the normalized velocity profile of the Newtonian fiuid does not depend on either the 
pressure gradient AP or the temperature T. 
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(b) local viscosity 

FIG. 5: Spatial variation of the local strain rates [for (a)] 
and the local melt viscosities [for (b)]. The profiles in the 
lower-half space are plotted. The profiles in the upper-half 
space are antisymmetric for the strain rate and symmetric for 
the viscosity. 



polymer melt its dependence is quite non-linear. Espe- 
cially, at a low temperature T = 0.2, Vq increases more 
than tenfold while the pressure gradient increases only 
threefold; = 0.024 to 0.31 for AP = 1 x lO"'^ to 
3 X lO^'^. 



In Fig. [5l we show the spatial variations of the local 
strain rates and local viscosities of the melts. Here, the 
local strain rates are time-averaged as in the velocity pro- 
files. The local viscosities are calculated by dividing the 
time averages of the local stresses by those of the local 
strain rates. It is seen that the shear thinning occurs 
near the plate; i.e., the local viscosities becomes thin- 
ner as the local strain rates near the plate increase. The 
shear thinning is enhanced at a low temperature and a 
high pressure gradient. For a melt with a low tempera- 
ture T=0.2 and a high pressure gradient AP=5xlO~^, 
local viscosity near the plate decreases to less than 3% of 
that in the middle. Thus, the velocity gradient becomes 
quite large compared to that for the middle between the 
plates, so that the flatter velocity profile shown in Fig. [3] 
is produced. 

The microscopic configurations of polymer chains are 
also investigated using the bead distribution function 175 
and structure factor S defined below, which are calcu- 
lated from the microscopic configurations of bead par- 
ticles obtained at each MD cell. The bead distribution 
function gs(r) and structure factor S{q) are defined by 

with R% = N\,-^ Y.^2i Rt and 

where k is the index to represent the fcth polymer chain 
and i2f represents the position of the ith bead parti- 
cle on the fcth polymer chain in each MD cell. Here, 
the assembled average {■ ■ ■)k is calculated by the average 
over all the configurations of polymer chains obtained 
at i = 100000 -I- lOOOn {n = 0, ■ • ■ , 100) in each MD 
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(b) y = 775 

FIG. 6: The local bead distribution function gs{r) defined 
by Eq. (|17p in the Vx ~ r-y plane [r^ = 0) and local structure 
factor S{q) defined by Eq. (I18|) in the qx — qy plane (g^ — 0) 
for the melt for Case B at y=25 [for (a)], j/=375 [for (b)] and 
y=775 [for (c)]. The contour lines show the values 0.01 + 
0.02fc for Qs and 0.05 + 0.2k for S with fc = 0, 1, • ■ ■ ,4 from 
outer to inner. 



cell. In Fig. [51 we plot the bead distribution function 
gsand structure factor S for the case with T = 0.2 and 
AP = 3 X 10^^ (Case B) at different positions. It is seen 
that the coherent structure is produced near the plate, 
while the incoherent structure is produced in the middle 
between the plates. The polymer chains are elongated 
in the flow direction (x-axis) near the plate, while in the 
middle between the plates bead particles are distributed 
rather randomly around the center of mass of each poly- 
mer chain. 



C. Oscillating plates 

The lower and upper plates begin to oscillate with a 
constant angular frequency luq a-t a time t = as, respec- 
tively, 

Vw = ±vo cos(woi), (19) 



TABLE III: The parameter values for the problem of oscillat- 
ing plates and the thickness of the boundary layer produced 
in each case (which is also shown in Figs. [5] and [6}. 





To 


27r/t<>o 


To 


H 


L 


Case I 


0.2 


1024 


0.5 


1600 


750 


Case II 


0.2 


256 


0.5 


800 


200 


Case III 


0.4 


256 


0.1 


800 


250 


Case IV 


0.2 


256 


0.1 


800 


550 




(a) Polymer Melt (b) Navier Stokes 



FIG. 7: The snapshots of velocity profiles in oscillating plates 
at LOot = 507r + 6, e/n^O, 1/2, 3/4, 1, 3/2, and 7/4, for 
lio ~ 27r/1024 and ro=0.5. (a) The result for the polymer 
melt with To — 0.2 [Case I], (b) The result for the Newtonian 
fiuid. The vertical axis shows the height y and the horizontal 
axis shows the velocity Vx- K represents the thickness of the 
boundary layer, in which the amplitude of the local oscillating 
velocity is more than 1 % of that of the oscillating plate, 
\vx\/vQ > 0.01. 



with vq = TqluqH. Here Fq represents the amplitude of 
strain of the system. Thus, the strain of the system F(t) 
in Eq. (fT2)) is written as 

r{t) ^Tosm{ujot). (20) 

We perform the hybrid simulations for the parameters 
listed in Table IIIII The distance between the plates 2H 
is divided into sixty-four slits. The mesh system is the 
same as that shown in Fig[TJb). Thus sixty- four MD cells 
are associated with the centers of each slit. The saving 
factor Ax/Imb is 5 for Case I in Table Hill and 2.5 for 
Case II-IV in Table HIl 

Figures [71 and [H show the instantaneous velocity pro- 
files between the oscillating plates. The thickness of the 
viscous boundary layer near the oscillating plate l^, which 
is defined by the specification that the amplitude of the 
local oscillating velocity is 1% of that of the oscillating 
plate, i.e., \vx{y = ^i/)|/wo — 0.01, is also shown in each 
figure. Figure [7l shows the comparison of the melt for 
Case I and the Newtonian fluid with a constant viscosity 
i> = 228, which corresponds to the dynamics viscosity of 
the model polymer melt at ujq = 27r/1024 for small strain 
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(a) ro=0.5, T=0.2 



(b) ro=0.1, f =0.4 



(c) ro=0.1, T=0.2 



FIG. 8: The snapshots of velocity profiles in oscillating plates at LJot — IOOtt + 9, 9/n=0, 1/2, 3/4, 1, 3/2, and 7/4, for 
cDo = 27r/256. (a) The result for Case II, (b) that for Case III, and (c) that for Case IV. See also the caption for Fig. [S] 



rates. The dynamics viscosity is calculated by G" /ujq 
(where G" is the loss modulus in the linear regime). The 
velocity profiles are quite different from each other. It is 
seen that the boundary layer of the melt is much thinner 
than that of the Newtonian fluid. (The thickness of the 
boundary layer may be expressed as l^, oc a/ Pq /wq for 
the Newtonian fluid.) This is caused by shear thinning; 
the local loss modulus near the boundary is, as we see 
in Fig. [ini much smaller than that of the linear regime. 
Figure [5] shows the comparison of the velocity profiles for 
Cases II-IV in Table IIIII The comparison of (a) and (c) 
shows the effect of the amplitude of the oscillating veloc- 
ity and that of (b) and (c) shows the effect of the melt 
temperature. The viscosity becomes smaller as the tem- 
perature or strain rate increase, and thus the boundary 
layers become thinner in Case II and Case III compared 
with Case IV. 

We also measure the "local" visco-elastic properties in 
terms of the "local" storage modulus G' and loss mod- 
ulus G". It should be noted that the local macroscopic 
quantities oscillate with different phase retardations at 
each separate point. The local moduli are calculated as 
follows: The discrete Fourier transforms of the tempo- 
ral evolutions of the strain 7, ^{t,y) = /g 'y{t',y)dt', and 
shear stress a^y during the steady oscillation states are 
performed, and are written as 

1 ^ 

ffl-^^E^ie-''^'"-')"^-^^/^ (fc = l,...,iV), (21) 

with gl^=g{nAt,lAx) {n~l,...,N and Z=0,...,64) , where 
g represents the strain or shear stress (e.g., g—^ or (Jxy)- 
Hereafter the subscript k represents the mode in Fourier 
space. The time evolution of the local strain at y=y^ 
may be expressed by using the Fourier coefficients for 
the mode of oscillation of the plate ko, fco=l-l-(wo/27r)iV, 
as 

V(i) = |7l'cos(c^oi + '5'), (22) 



with |7|' = ^Re(27[o))2+Im(27[o)2 ^nd = 

tan~^(Im(7[Q)/Re(7[g)). The time evolution of the lo- 
cal shear stress for the mode ko can also be expressed 
as 

<^iy(t) = cos{ujot + (5') - sm{uJot + S'), (23) 
where 

a[ = Re{2ai^ ) cos 5' -I- Im(24o ) sin , (24a) 
cr^ = Im(24 J cos 5' -Re(24„)sin(5'. (24b) 

Thus, the local storage modulus G' and loss modulus 
G" are obtained, respectively, as G"(?/')=(tJ /I7I' and 
{y'')=o'2/\j\'' ■ We note that, in the Fourier trans- 
formations, the signals for Sloq Sind 5luq are also de- 
tected, although their contributions are much smaller 
than that for ajQ. In Table IIV| we show the ampli- 
tudes of the harmonic contributions at the frequency Sluq 
and 5wo for the local stresses. |cr^j^(tt')| is the spector of 
local stress in the Fourier space, which is obtained by 
Wiy{Lu)\ = ^a[{ky+a(,{k)^ with fc = 1 + {uj/2tt)N. 
Here, a[ 2{k) are calculated via Eq. ([M)) by replacing 
fco with k. Table IIVI represents the amplitude of non- 
linear response with respect to the frequency. It is seen 
that the non-linear response is enhanced in the boundary 
layer near the oscillating plate. 

Figure [5] shows the spatial variations of the amplitude 
of the local strain |7'| and local phase retardation 6^ in 
Eq. for Case I and Case II. The profiles of the am- 
plitude of the local strain I7I between the plates are sim- 
ilar in Case I and II. These amplitude profiles increase 
rapidly in the boundary layer near the oscillating plate, 
and have values much greater than unity, I7I ~ 10, in the 
vicinity of the plate. Thus, the polymer chains are quite 
deformed near the oscillating plates, although the ampli- 
tude of strain to which the system is subjected via the 
oscillating plates is not so high; Fq = 0.5 in both cases. 

Figure [TO] shows the spatial variations of the local stor- 
age modulus and loss modulus for Case I and Case II. 
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TABLE IV: The harmonic contributions of spectors of the 
local stress |CT^j^(ti;)| for frequency Sluq and 5iUo in Case I and 
II. 



Case I 



y 




!o-^a(3i^o)| 
\(j^y{ujo)\ 


\(Jxy(bljJo)\ 
\c^xy{^o)\ 


25 


0.938 


14% 


4% 


125 


0.428 


15% 


6% 


525 


0.106 


10% 


3% 


725 


0.075 


7% 


2% 




U.UDD 


CO/ 

o/o 
Case II 


U.z/o 


y 




|o-^h(3i^o)| 
\a^y{ujo)\ 


\axy{^i^o)\ 

\o'xy{i^o)\ 


12.5 


2.097 


10% 


3% 


62.5 


0.594 


13% 


5% 


112.5 


0.294 


13% 


3% 


212.5 


0.139 


9% 


1% 


362.5 


0.078 


3% 


1% 




ItI (■) l7l (■) 

(a) Case I (b) Case II 

FIG. 9: The spatial variations of I7I and 5 in Eq. ([22|l for 
(a) Case I and (b) Case II. The black square ■ shows I7I and 
the diamond O shows J/tt. 



Shear thinning is seen near the plate; both moduh G' 
and G" decrease near the oscillating plate. In the close 
vicinity of the oscillating plate, the storage modulus G' 
is much smaller than the loss modulus G", G" ^ G" . 
Hence, the melt behaves as a viscous fluid. The stor- 
age modulus grows rapidly with the distance from the 
oscillating plate, and the visco-elastic crossover occurs 
at y 1000 for Case I and at y ^ 200 for Case II. 
Both moduli attain their linear values, which are shown 
as dashed and dot-dashed lines in the figures, for a dis- 
tance that is far from the oscillating plate where the local 
strains are less than about 2 % (See also Fig. [Q]). Thus, 
the local rheology of the melt can be divided into three 
regimes, i.e., the viscous fluid, visco-elastic liquid, and 
visco-elastic solid regimes. 




G", G" 
(b) Case II 

FIG. 10; The spatial variations of the local moduli G" (■) 
and G" (O) for (a) Case I and (b) Case II. The dashed and 
dash-dotted lines show the values of & and G" for the linear 
regime (G" = 1.7 and G" = 1.4 for Case I and G' = 3.1 
and G" = 1.3 for Case II), respectively. The linear moduli 
are calculated by the non-equilibrium MD simulations with 
small strains 0.005 < I7I < 0.01. The left arrows on the right- 
side vertical axis show the positions where the local Deborah 
numbers, shown in Fig. 1101 are equal to unity and the position 
where the local strain is I7I = 2%. 



These regimes may be also characterized by the two 
"local" Deborah numbers. One is defined by the local 
Rouse relaxation time of the melt and the angular 
frequency of the plate t^Oj De^—ujQTfj, and the other is 
defined by the local a relaxation time Tq, and the angular 
frequency ujq, De°'=aJoTa. Note that the local Rouse and 
a relaxation times vary according to the local strain rate 
7, r = r(7). Figure [TT] shows the spatial variation of 
the local Deborah number De^ and De" , where the local 
relaxation times and are calculated by substituting 
the values of I7I', which are obtained by Eq. pTjl and 
equation below Eq. (f^^ . into the fitting functions for 
the relaxation times for the simple shear flows obtained 
in Ref. Q- In Fig. [TOl the positions at which the local 
Deborah numbers become equal to unity are shown by 
the left arrows. It is seen that the melt behaves as a 
viscous fluid, G" ^ G', for De^ < 1, while the visco- 
elastic properties become pronounced for De^ > 1. This 
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FIG. 11: The spatial variations of the local Deborah numbers 
(defined via the Rouse relaxation time tr and the a relaxation 
time t") L>e^ = ujotr (Q) and _De" = uor" (■). (a) The 
result for Case I and (b) that for Case II. The dashed line 
represents De = 1. 



feature is also consistent with the rheology diagram for a 
model polymer melt obtained from Ref. [23. It is also seen 
that the crossovers from the liquid-like regime, G" > G', 
to the solid- like regime, G" > G", take place at De°' ~ 1. 



IV. SUMMARY AND OUTLOOK 

The behaviors of supercooled polymer melt in creep 
and recovery, pressure-driven flow and oscillating flows 
between the parallel plates are investigated numerically 
by using a hybrid simulation of MD and CFD. In the 
present hybrid method, the memories of molecular con- 
figurations of local fluid elements are traced at the micro- 
scopic level so that the visco-elastic motion of the melt 
is correctly reproduced. 

The flow profiles of the melts are quite different from 
those of the Newtonian fluid. In the creep simulation, wc 
demonstrate the simple visco-elastic motion of the melt. 
The nonlinear time evolution of strain of the system at 
the beginning of creep motion, such as the delayed elastic 
deformation, is reproduced. After removing the stresses 
on the plates, the polymer melt recovers by an elastic 
strain stored in the creep motion. We also compare the 
result given by the present hybrid method with that given 
by the Jeffreys model for a constitutive relation. It is seen 
that the result from the hybrid method is fitted well with 
that given by the Jeffreys model, although the fitting 
parameters in the Jeffreys model differ from each other 
in the creep and recovery motions since the mechanical 
property of the melt is quite sensitive to motion states. 
(See Fig. [21) In pressure-driven flow, shear thinning oc- 
curs near the plates, where the local shear rates are much 
larger than those in the middle of the plates. Shear thin- 
ning is enhanced as the pressure gradient increases or 
the temperature decreases, so that the velocity profile is 
increasingly flatter rather than parabolic at the middle 
between the plates. (See Figs. [Hand [51) We also inves- 



tigate the microscopic configurations of polymer chains 
at local points. The polymer chains are quite elongated 
in the flow direction near the plate, while in the middle 
between the plates, conflgurations of polymer chains are 
rather incoherent. (See Fig. [6l) For oscillating plates, 
we clarify that the viscosity of the melt becomes thin 
near the plates, and thus the boundary layer of the melt 
also becomes much thinner than that of the Newtonian 
fluid. (See Figs. [7land[8l) The local rheological prop- 
erties of the melt also vary considerably in the viscous 
boundary layer, so that three different regimes form be- 
tween the oscillating plates, i.e., the viscous fluid, visco- 
elastic liquid, and visco-elastic solid regimes. It is also 
found that, in the viscous fluid regime, the local Debo- 
rah number defined via the Rouse relaxation time and 
the angular frequency of the plate is approximately less 
than unity, ZJe^ < 1. The crossover between the liquid- 
like and solid-like regimes takes place around the position 
where the local Deborah number (defined via the a relax- 
ation time and the angular frequency) is approximately 
equal to unity, De" - 1. (See Figs. [iniand[ni) 

In the present hybrid method, the long-range correla- 
tions of the macroscopic quantities that are difficult to 
treat at the MD level are involved at the CFD level via 
the macroscopic momentum transport equation. Thus, 
although each MD simulation is performed independently 
at each time step of CFD, the polymers in different MD 
cells are also correlated with each other via the macro- 
scopic momentum transports. The complicated mechan- 
ical properties depending on the microscopic configura- 
tions of polymer chains for which it is difficult to con- 
struct the model constitutive relation are also correctly 
involved in the CFD simulation, since the microscopic 
motions of polymers are resolved in the MD simulations 
associated with each mesh node of the CFD calculation 
according to the local flow fleld. Hence, the present hy- 
brid method is expected to be applicable both to prob- 
lems for the large-scale flows which are out of range of 
the MD simulation, and the complex flows, which do not 
have any known model constitutive relations. 

Compared with running a full MD simulation, the 
present hybrid simulation can save computation time as 
to the spatial domain, although it does not accelerate 
computation time with respect to the the temporal do- 
main. The efficiency of this hybrid simulation is repre- 
sented by a saving factor defined by the ratio of the mesh 
size of the CFD simulation Ax to the cell size of the MD 
simulation Zmd, Ax/Iub- In the present simulations, the 
saving factors are Aa;/ZMD=4.7 for the creep and recov- 
ery, Aa;/ZMD=5 for pressure-driven flow and Case I in the 
oscillation problem, and Ax/?md=2.5 for Cases II-IV in 
the oscillation problem. In addition to the saving factor 
Ax/Imd, the present hybrid method has also unique in 
that it is quite suitable as a parallel computational al- 
gorithm since the MD simulations associated with each 
mesh of the CFD, which represent a large part of the 
total simulation time, are performed independently. 

We carry out benchmarks for parallel computations for 
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FIG. 12: Benchmarks for parallel computations using the 
hybrid simulation method. A'^cpu is the number of CPUs and 
tNcpxj the computational time of parallel computation with 
A^CPU CPUs. The speed of parallel computation with A^cpu 
CPUs is measured by the ratio of the computational time 
using one CPU to that with A'^cpu CPUs. Parallel computa- 
tions for A^cpu=4, 16, 32, 64, 128, and 256 were performed on 
Fujitsu HX600 at T2K Open Supercomputer at Kyoto Uni- 
versity. 

the hybrid simulation method. In the benchmarks, the 
creep simulations with 2H — 12800 and iro—Q.l are per- 
formed for 1,000 of CFD, i = [0,1000]. The distance 
between the plates 2H is divided into 257 shts. Thus 
256 MD cells with ^md = 10 are associated with each 
node of the shts. The parallehzation algorithm is ap- 
phed to the process to calculate local stresses by MD 
simulation. Hence, in the parallel computations with 
-^CPU CPUs, each CPU is assigned to the calculations 
of (256/A'cpu) local stresses at each time step of the 
CFD. Figure [T^] shows the result of the benchmarks. 



The parallel computations are performed nearly ideal 
until the number of CPU TVcpuis one-hundred twenty- 
eight, A^cpu = 128. The efficiency of parallel computa- 
tion e defined by e = {ti/tjs[^p-fj)/Ncpi!i where tj^^^^ is 
the computational time for iVcpu parallehzation, is more 
than 99% up to Nqpu = 128 and 90% for Nqpu = 256. 
Parallehzation efficiency decreases a bit at iVcpu=256 in 
the present benchmark since CFD calculation and com- 
munication between CPUs increases relative to the total 
computation. However, for the simulations of large-scale 
and slow dynamics, we can improve the parallehzation 
efhciency even for a large number of CPUs by setting 
the slit size and time-step size of the CFD calculation. 
Ax and At, large while keeping the ratios Ax/Imb and 
At/tMB constant (because the ratio of CFD calculation 
computation time to the total time decreases.) Although 
the total computation time dose indeed increases, the 
efficiency of the hybrid simulation compared to the full 
MD simulation does not change. Thus we can expect to 
carry out quite high-performance parallel computations 
even for thousands or tens of thousands of CPUs. 

In the present paper we are involved with only one- 
dimensional problems. Future work includes the develop- 
ment of the hybrid method for two- or three-dimensional 
simulations. In order to extend the present method to 
the two- or three-dimensional problems, we must treat 
the advection of memory of polymer chain configurations 
on the local fluid element. For two-dimensional or three- 
dimensional flows of visco-elastic fluids, many novel sim- 
ulation methods have been proposed,2i^i2Ii22ii22 Incorpo- 
rating the present hybrid technique into well-established 
simulation methods for the visco-elastic flows represents 
an important direction for future research. The exten- 
sion of this technique to the coupling of macroscopic heat 
and mass transfers is also important to the investigation 
of various physical problems involved in polymer melts. 
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